Extract bubble from simulations

In [1]:
import numpy as np
import scipy as scp
from scipy import optimize
import matplotlib.pyplot as plt

Simulation Parameters

In [2]:
nLat = 4096
nSims = 100
minSim = 0
nCols = 4

phi_initial = np.pi
nu = 2.*10**(-3)
lamb = 1.5; print('lamb = ', lamb)
m2eff = 4. * nu * (- 1. + lamb**2); print('m2eff = ', m2eff)
lenLat = 2 * 50. / np.sqrt(2. * nu); print('lenLat = ', lenLat)
phi0 = 2. * np.pi / 4.; print('phi0 = ', phi0)

lamb =  1.5
m2eff =  0.01
lenLat =  1581.1388300841897
phi0 =  1.5707963267948966
In [3]:
pickle_location = '/gpfs/dpirvu/thick_wall_average_bubble/'
suffix = '_for_phi0{:.4f}'.format(phi0)+'_lambda'+str(lamb)+'_len{:.4f}'.format(lenLat)+'_x'+str(nLat)

def bubbles_file(min, max):
    return pickle_location+'bubbles_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def sim_location(sim):
    return '/gpfs/dpirvu/sims/x'+str(nLat)+'_phi0{:.4f}'.format(phi0)+'_lambda{:.4f}'.format(lamb)+'_sim'+str(sim)+'_fields.dat'

def V(phi):
    return ( -np.cos(phi) + 0.5 * lamb**2 * np.sin(phi)**2 ) * 4 * nu
def F(x):
    return V(x) - V(phi_initial)

phi_upper_bound = optimize.fsolve(F, 5)[0]; print(phi_upper_bound)
phi_lower_bound = optimize.fsolve(F, 0.5)[0]; print(phi_lower_bound)

4.823729994725654
1.4594553124539482

Extract data from files

In [4]:
def extract_data(filename, col):
#    print(filename)
    infile = open(filename,'r')
    lines = infile.readlines()
    field_values = [float(line.split()[col]) for line in lines[3:]]
    infile.close()
    return field_values

def check_decay(simulation):
    right_phi = sum([1 for x in simulation[-1] if x > phi_upper_bound])
    left_phi = sum([1 for x in simulation[-1] if x < phi_lower_bound])
    if right_phi > nLat*0.1 and left_phi < nLat*0.02: return 0
    elif left_phi > nLat*0.1 and right_phi < nLat*0.02: return 1
    else: return 2

def time_at_fraction(bubble, frac, limit):
    T, N = len(bubble), len(bubble[0])
    right_phi_x = [np.sum([1 for x in slice if x >= limit]) for slice in bubble]
    time_list = [t if (right_phi_x[t] <= N*frac) else 0 for t in range(T)]
    return next((t for t in time_list[::-1] if t != 0), 0)

def triage(sim):
    y = extract_data(sim_location(sim), 0)
    nT = len(y)//nLat
    tmin = 3000
    if nT > tmin:
        outcome = check_decay(np.reshape(y, (nT, nLat)))
        if outcome != 2:
            data = [extract_data(sim_location(sim), col) for col in range(nCols)]
            try:
                data = [np.reshape(data[col], (nT, nLat)) for col in range(nCols)]
            except ValueError:
                return None
            if outcome == 1:
                data[0] = [2*phi_initial - i for i in data[0]]
                if nCols > 1:
                    data[1] = [-j for j in data[1]]
            tdecap = time_at_fraction(data[0], 0.01, phi_upper_bound)
            tdecap0 = tdecap + 2000
            data = [data[i][max(tdecap-tmin, 0):tdecap0] for i in range(len(data))]
            print('sim', sim, ', duration ', nT)
            return data, sim

def alldata():
    all_rsp_data, sims_to_keep = [], []
    for sim in range(minSim, nSims):
        if sim == nSims // 2:
            print('Halfway through.')
        a = triage(sim)
        if a is not None:
            data, sim = a
            all_rsp_data.append(data)
            sims_to_keep.append(sim)
    return all_rsp_data, sims_to_keep

def checkdata():
    all_rsp_data = []
    for sim in range(minSim, nSims):
        if sim == nSims // 2:
            print('Halfway through.')
        y = extract_data(sim_location(sim), 0)
        nT = len(y)//nLat
        data = [np.reshape(extract_data(sim_location(sim), col), (nT, nLat)) for col in range(nCols)]
        all_rsp_data.append(data)
    return all_rsp_data
In [5]:
all_data, sims_to_keep = alldata()
#all_data = checkdata()
print(np.shape(all_data))

sim 2 , duration  3074
sim 5 , duration  10211
sim 11 , duration  5094
sim 12 , duration  6026
sim 15 , duration  6176
sim 17 , duration  4763
sim 19 , duration  14207
sim 24 , duration  21370
sim 25 , duration  4272
sim 28 , duration  15037
sim 29 , duration  6793
sim 31 , duration  4428
sim 33 , duration  18701
sim 34 , duration  15176
sim 36 , duration  13589
sim 46 , duration  3957
sim 47 , duration  17181
Halfway through.
sim 52 , duration  17022
sim 55 , duration  12356
sim 56 , duration  10226
sim 59 , duration  9004
sim 67 , duration  17537
sim 70 , duration  3257
sim 72 , duration  3782
sim 74 , duration  4340
sim 77 , duration  5483
sim 79 , duration  18717
sim 81 , duration  18742
sim 90 , duration  12326
sim 94 , duration  10865
sim 96 , duration  6348
sim 97 , duration  15265
sim 98 , duration  11031
sim 99 , duration  5753
(34, 4)
In [6]:
np.save(bubbles_file(minSim, nSims), [all_data, sims_to_keep])
In [7]:
def plot_real_space_data(sim, col):
    fig, ax0 = plt.subplots(1, 1, figsize = (7, 5))
    im0 = ax0.imshow(all_data[sim][col], aspect='auto', interpolation='none', origin='lower')
    clb = plt.colorbar(im0, ax = ax0); clb.set_label(r'$\phi(x)$', labelpad=-48, y=1.08, rotation=0)
    ax0.set(xlabel = r'$x$', ylabel = r'$t$')#; ax1.set(xlabel = r'$x$', ylabel = r'$t$')
    plt.show()
    return

def plot_real_space_slice(sim, col, timeslice):
    slice = all_data[sim][col][timeslice]
    plt.figure(figsize = (50, 3))
    plt.plot(np.arange(len(slice)), slice, label=timeslice)
    labelLines(plt.gca().get_lines(), xvals=(0, nLat//2), align=False)
    plt.xlabel(r'$x$'); plt.ylabel(r'$\phi(x)$'); plt.legend(); plt.show()
    return
In [8]:
for sim in range(len(all_data))[::]:
    plot_real_space_data(sim, 0)